SARIMA

Ejercicio de Simulación

Vamos a hacer un ejercicio de simulación para ver como se identifica la componente estacional.

library(urca)
library(forecast)
library(tseries)
library(lmtest)
library(uroot)
library(fUnitRoots)

Loading required package: timeDate
Loading required package: timeSeries

Attaching package: ‘timeSeries’

The following object is masked from ‘package:zoo’:

    time<-

Loading required package: fBasics

Attaching package: ‘fUnitRoots’

The following objects are masked from ‘package:urca’:

    punitroot, qunitroot, unitrootTable

library(sarima)

Loading required package: stats4

Attaching package: ‘sarima’

The following object is masked from ‘package:stats’:

    spectrum

require("PolynomF")

Loading required package: PolynomF

###Simulación de un proceso con raíz unitaria estacional
#x11()
x <- ts(sarima::sim_sarima(n=144, model = list(iorder=0, siorder=1, nseasons=12, sigma2 = 1),n.start=24),frequency = 12)
plot(x)

acf(x,lag.max = 36)

monthplot(x)

nsdiffs(x)####Decreta cuantas diferencias estacional a través de la aplicación de 

[1] 1

###Algunas pruebas de raíces unitarias estacionales.


###diferencia estacional
Dx=diff(x,lag=12,differences = 1)###lag:periodo s.
plot(Dx)

acf(Dx,lag.max = 36)

monthplot(Dx)

nsdiffs(Dx)

[1] 0

####Simulación de un SAR
#x11()
x1 <- ts(sim_sarima(n=144, model = list(ar=c(rep(0,11),0.8)),n.start=24),frequency=12)
plot(x1)

acf(x1,lag.max = 36)

monthplot(x1)

nsdiffs(x1)

[1] 1

ndiffs(x1)

[1] 1

p <- polynom(c(1,c(rep(0,11),-0.8)))
solve(p)

 [1] -1.0187693+0.0000000i -0.8822801-0.5093846i
 [3] -0.8822801+0.5093846i -0.5093846-0.8822801i
 [5] -0.5093846+0.8822801i  0.0000000-1.0187693i
 [7]  0.0000000+1.0187693i  0.5093846-0.8822801i
 [9]  0.5093846+0.8822801i  0.8822801-0.5093846i
[11]  0.8822801+0.5093846i  1.0187693+0.0000000i

abs(solve(p))

 [1] 1.018769 1.018769 1.018769 1.018769 1.018769 1.018769
 [7] 1.018769 1.018769 1.018769 1.018769 1.018769 1.018769

###Note lo cerca que están la raíces de la no estacionariedad del proceso, por eso
####aunque si bien el proceso es estacionario, notamos hay una cercanía a 
####e tener una compoenete estacional.
####El anterior modelo puede escribirse como:
x2 <- ts(sim_sarima(n=144, model=list(sar=0.8, iorder=0, siorder=0, nseasons=12),n.start=24),frequency = 12)
plot(x2)

acf(x2, lag.max=48)

monthplot(x2)

nsdiffs(x2)

[1] 1

Ejemplo Pasajeros

Vamos a ver como se hace el modelamiento completo de la serie de pasajeros.

Iniciaremos con la transformación Box-Cox y las pruebas de raíces Unitarias

[1] 1

Call:
ar(x = lAirpassengers)

Coefficients:
      1        2        3        4        5        6  
 1.0013  -0.0959   0.0329  -0.0252   0.0231   0.0005  
      7        8        9       10       11       12  
-0.0157  -0.0762   0.1071  -0.0236   0.0675   0.4536  
     13  
-0.4854  

Order selected 13  sigma^2 estimated as  0.01336

    Augmented Dickey-Fuller Test

data:  lAirpassengers
Dickey-Fuller = -2.147, Lag order = 13, p-value =
0.5152
alternative hypothesis: stationary

Warning in fUnitRoots::adfTest(lAirpassengers, lags = 12, type = "nc") :
  p-value greater than printed p-value

Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: 3.7872
  P VALUE:
    0.99 

Description:
 Thu Jun  9 11:01:18 2022 by user: 


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.103864 -0.023652 -0.001455  0.022160  0.126649 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
z.lag.1       0.005439   0.001436   3.787 0.000241 ***
z.diff.lag1  -0.198813   0.072080  -2.758 0.006737 ** 
z.diff.lag2  -0.273730   0.073309  -3.734 0.000292 ***
z.diff.lag3  -0.233607   0.072258  -3.233 0.001589 ** 
z.diff.lag4  -0.293133   0.073926  -3.965 0.000126 ***
z.diff.lag5  -0.206562   0.072058  -2.867 0.004915 ** 
z.diff.lag6  -0.266919   0.071493  -3.734 0.000292 ***
z.diff.lag7  -0.234526   0.071847  -3.264 0.001436 ** 
z.diff.lag8  -0.327393   0.073197  -4.473 1.79e-05 ***
z.diff.lag9  -0.198455   0.073623  -2.696 0.008054 ** 
z.diff.lag10 -0.279931   0.072710  -3.850 0.000192 ***
z.diff.lag11 -0.176122   0.073011  -2.412 0.017394 *  
z.diff.lag12  0.627402   0.072683   8.632 3.34e-14 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04308 on 118 degrees of freedom
Multiple R-squared:  0.8564,    Adjusted R-squared:  0.8406 
F-statistic: 54.14 on 13 and 118 DF,  p-value: < 2.2e-16


Value of test-statistic is: 3.7872 

Critical values for test statistics: 
      1pct  5pct 10pct
tau1 -2.58 -1.95 -1.62


Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.5325
  P VALUE:
    0.7711 

Description:
 Thu Jun  9 11:01:18 2022 by user: 

Call:
ar(x = dlAirpassengers)

Coefficients:
      1        2        3        4        5        6  
-0.0993  -0.1226  -0.2824  -0.2875  -0.1184  -0.2406  
      7        8        9       10       11       12  
-0.1829  -0.2971  -0.1040  -0.2627  -0.1259   0.5734  
     13       14       15  
 0.0179  -0.1667   0.1200  

Order selected 15  sigma^2 estimated as  0.002988

Warning in tseries::adf.test(dlAirpassengers, k = 2) :
  p-value smaller than printed p-value

    Augmented Dickey-Fuller Test

data:  dlAirpassengers
Dickey-Fuller = -7.6877, Lag order = 2, p-value =
0.01
alternative hypothesis: stationary


Title:
 Augmented Dickey-Fuller Test

Test Results:
  PARAMETER:
    Lag Order: 12
  STATISTIC:
    Dickey-Fuller: -1.2228
  P VALUE:
    0.2243 

Description:
 Thu Jun  9 11:01:18 2022 by user: 

Identificación de la componente ARMA estacional y la componente ARMA ordinaria

####################################
library(uroot)
require(forecast)
######Diferencia Estacional(continuación AirPassengers)#######
monthplot(dlAirpassengers)
nsdiffs(dlAirpassengers)

[1] 1

nsdiffs(AirPassengers)

[1] 1

DdlAirpassengers=diff(dlAirpassengers,lag=12)###lag=s
#x11()
par(mfrow=c(2,1))

plot(dlAirpassengers)
plot(DdlAirpassengers)

monthplot(DdlAirpassengers)
nsdiffs(DdlAirpassengers)

[1] 0

##Autocorrelogramas
#x11()
acf(DdlAirpassengers)

acf(DdlAirpassengers,lag.max = 48, ci.type='ma')# q=0,1, Q=0,1
pacf(DdlAirpassengers,lag.max = 48) # p=0,1,2,...,9, P=0,1

#SARIMA(p=0,d=1,q=1)x(P=0,D=1,Q=1)s=12

Ajuste del Modelo y Análisis de Residuales

##Ajuste del modelo
###Arima Estacional o SARIMA(p=0,d=1,q=3)x(P=0,D=1,Q=1)s=12 con transformación logaritmica

#Modelo MA(1) estacional
modelo = Arima(AirPassengers, c(0, 1, 1),seasonal = list(order = c(0, 1, 1), period = 12),lambda = 0)
coeftest(modelo)

z test of coefficients:

      Estimate Std. Error z value  Pr(>|z|)    
ma1  -0.401828   0.089644 -4.4825 7.378e-06 ***
sma1 -0.556945   0.073100 -7.6190 2.557e-14 ***
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

modeloalter= Arima(AirPassengers, c(1, 1, 0),seasonal = list(order = c(1, 1, 0), period = 12),lambda = 0)

## Análisis de residuales
#x11()
residuales <- modelo$residuals
plot(residuales)

acf(residuales,lag.max = 24)

pacf(residuales,lag.max = 24)
#Test de autocorrelaci?n
Box.test(residuales, lag = (length(residuales)/4), type = "Ljung-Box", fitdf = 2)

    Box-Ljung test

data:  residuales
X-squared = 37.874, df = 34, p-value = 0.2969

######Análisis de Outliers
#Test de normalidad
jarque.bera.test(residuales)

    Jarque Bera Test

data:  residuales
X-squared = 5.2265, df = 2, p-value = 0.0733

###Estad?ticas CUSUM
res=residuales
cum=cumsum(res)/sd(res)
N=length(res)
cumq=cumsum(res^2)/sum(res^2)
Af=0.948 ###Cuantil del 95% para la estad?stica cusum
co=0.14422####Valor del cuantil aproximado para cusumsq para n/2
LS=Af*sqrt(N)+2*Af*c(1:length(res))/sqrt(N)
LI=-LS
LQS=co+(1:length(res))/N
LQI=-co+(1:length(res))/N
par(mfrow=c(2,1))

plot(cum,type="l",ylim=c(min(LI),max(LS)),xlab="t",ylab="",main="CUSUM")
lines(LS,type="S",col="red")

lines(LI,type="S",col="red")
#CUSUM Square
plot(cumq,type="l",xlab="t",ylab="",main="CUSUMSQ")                      

lines(LQS,type="S",col="red")                                                                           
lines(LQI,type="S",col="red")

Pronóstico

ecm

[1] 342.0955

Vale la pena decir que si los residuales son no correlacionados pero ellos no tiene distribución normal, se puede usar los intervalos bootatrap ver sección 5.5 libro fpp3. Esto se hace añadiendo bootstrap=TRUE en la función forecast().